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Abstract — Optimal power flow (OPF) is an important problem 
for power generation and it is in general non-convex. With 
the employment of renewable energy, it will be desirable if 
OPF can be solved very efficiently so its solution can be used 
in real time. With some special network structure, e.g. trees, 
the problem has been shown to have a zero duality gap and 
the convex dual problem yields the optimal solution. In this 
paper, we propose a primal and a dual algorithm to coordinate 
the smaller subproblems decomposed from the convexifled OPF. 
We can arrange the subproblems to be solved sequentially and 
cumulatively in a central node or solved in parallel in distributed 
nodes. We test the algorithms on IEEE radial distribution test 
feeders, some random tree-structured networks, and the IEEE 
transmission system benchmarks. Simulation results show that 
the computation time can be improved dramatically with our 
algorithms over the centralized approach of solving the problem 
without decomposition, especially in tree-structured problems. 
The computation time grows linearly with the problem size with 
the cumulative approach while the distributed one can have size- 
independent computation time. 

I. Introduction 

AN ELECTRIC power system is the main facility to 
distribute electricity in modem societies. It is a network 
connecting power supplies (e.g., thermoelectric generators and 
turbine steam engines) to consumers. A power grid is generally 
composed of several subsystems: generation, transmission, 
substation, distribution, and consumers. The generators gen- 
erate power which is delivered to the substations at high 
voltages through the transmission network. The power voltage 
is stepped down and then distributed to the consumers via 
the distribution networks. In a typical power system, a few 
hundreds of generators interconnect to several hundreds of 
substations. The substations distribute the power to millions of 
consumers with relatively simpler radial networks with tree- 
like structures. 

In the past, research on power systems mainly focused 
on the core of the network, i.e., from the generation, via 
transmission, to the substations. All of the control, planning 
and optimization was done by a single entity (e.g. an ISO). 
With the integration of renewal energy and energy storage, 
self-healing ability, and demand response, the focus is shifted 
toward the consumer side, i.e. distribution networks, and this 
new paradigm is called the smart grid [T|. 

The optimal power flow (OPF) is one of the most important 
problems in power engineering and it aims to minimize the 
generation cost subject to demand constraints and the network 
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physical constraints, e.g. bus voltage limits, bus power limits, 
thermal line constraint, etc. Due to the quadratic relations 
between voltage and power, OPF is non-convex. In general, 
heuristic approaches have been employed to solve the OPF 
but they are not guaranteed to yield the optimal solution. To 
simplify the calculation, with assumptions on lossless power 
line, constant voltage and small voltage angles, OPF can 
be linearized and this approximation is also called DC-OPF, 
which is not accurate under all circumstances Q. For the 
complex OPF, |3| suggested solving the problem in its dual 
form and studied the conditions of the power network with 
zero duality gap. In Q, it was shown that the duality gap 
is always zero for network structures such as trees which 
model distribution networks well. |5|, as an independent work 
of this paper, decomposes the OPF in terms of cycles and 
branches and formulates the problem as an second-order cone 
program for tree networks which is equivalent to that given in 
In traditional power systems, OPF is mainly for planning 
purpose. For example, it is used to determine the system state 
in the day-ahead market with the given system information. 
In the smart grid paradigm, due to highly intermittent nature 
of the renewable, the later the prediction is made, the more 
reliable it is. If OPF can be solved very efficiently, we may 
solve the OPF in real time thus mitigating some of the 
unpredictability. 

We aim at solving OPF efficiently. When the system size 
(e.g. the number of buses) increases, solving the problem 
in a centralized manner is not practical (this will be veri- 
fied in the simulation). One possible way is to tackle the 
problem distributedly by coordinating several entities in the 
system, each of which handle part of the problem and their 
collaborative effort solves the whole problem. To do this, a 
communication protocol is needed to define what information 
should be conveyed among the entities. We can learn from 
the networking protocol development to design a commu- 
nication protocol for OPF. The earliest form of protocols 
for the Internet was proposed in 70's. They were designed 
to handle the increasing volume of traffic sent over the 
Internet in an ad hoc manner In 1998, Kelly et al. studied 
Transmission Control Protocol (TCP), which is one of the 
core protocols in TCP/IP |7|. They showed that TCP can be 
analyzed with a fundamental optimization problem for rate 
control and the algorithms developed from the optimization 
fit the ad hoc designed variants of TCP. This lays down a 
framework to design communication protocols for complex 
systems with reasoning. In this framework, we start with an 
optimization problem representing the system. By optimization 
decomposition |j8], the problem is decomposed into (sim- 
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pier) subproblems which can be solved by different entities 
in the system independently. The coordination between the 
subproblems define the communication protocols (i.e. what 
and how the data exchange between the entities). ||9| shows 
that many problems in communications and networking can 
be cast under this framework and protocols can be designed 
through primal and dual decomposition. In this paper, we study 
OPF by decomposing it into subproblems with primal and 
dual decomposition. Then we propose the primal and dual 
algorithms, respectively, to solve OPF in a distributed manner 
and the algorithms determine the communication protocols. 
Our algorithms do not assume the existence of a communica- 
tion overlay with topology different from the power network. 
In other words, a bus only needs to communicate with its 
one-hop neighbors in the power network. The algorithms can 
employed to any power network as long as the strong duality 
holds. We test the algorithms on IEEE radial distribution test 
feeders, some random tree-structured networks, and the IEEE 
transmission system benchmarks. Simulation results show that 
the algorithms can solve OPF distributively and they are 
scalable, especially when applied to the distribution network. 
If we apply the algorithm distributedly, the computational 
time is independent of the number of buses. If we apply the 
algorithms (with the problem decomposed) in a central node, 
the computational time grows linearly with the number of 
buses. 

The rest of this paper is organized as follows. In Section [H] 
we give the OPF formulation and the necessary background. 



Section III describes the primal and dual algorithms and the 
mechanism to recover the optimal voltage from the results of 
the algorithms. We illustrate the algorithms with two examples 



in Section M and present the simulation results in Section VI 



In Section [VII| we discuss the characteristics of the algorithms 



and conclude the paper in Section VIII 



II. Preliminaries 
A. Problem Formulation 

Assume that there are n buses in the power network. For 
buses i and k, i ^ k means that they are connected by a power 
line and i ^ k otherwise. Let zik and yik be the complex 
impedance and admittance between i and fc, respectively, and 
we have yn^ — j^. We denote Y = (F^fe, 1 < z, fc < n) as the 
admittance matrix, where 



Va if i^k 
-yik if i ~ fc 
if i 7^ fc. 



Let V = iVi,V2, . . . , KO^ G C" and i = (/i, /a, . . . , G 
C" be the voltage and current vectors, respectively. By Ohm's 
Law and Kirchoff's Current Law, we have i = Yv. The 
apparent power injected at bus i is Si — Pi + jQi — Vilf^ , 
where Pi and Qi are the real and reactive power, respectively, 
and H means Hermitian transpose. We have the real power 
vector p = {Pi,P2, . . . ,Pn)'^ = Re{diag(vv^Y^)}, where 
diag(vv^Y^) forms a diagonal matrix whose diagonal is 
We define the cost function of Bus i as costi{Pi) = 



Ci2P'l + CiiPi + Cio, where Cio,Qi,Ci2 S K and Ci2 > 0,Vi. 
OPF can be stated as 



minimize cost^ {Pi ] 



subject to 



]^<\v,\<v„yi 
p^<p. <P;-yi 

Ptk < Ptk^y-i, k 

p = Re{diag(vv^Y^)} 



(la) 

(lb) 
(Ic) 
(Id) 
(le) 



where Vi, Vi, Pi, Pi, and Pik are the lower and upper voltage 
limits of bus i, the lower and upper power limits of bus i, and 
the real power flow Umit between buses i and fc, respectively. 
Eq. (lb I is the nodal voltage constraint limiting the magnitude 



of bus voltage. Eq. ( Ic i is the nodal power constraint limiting 



the real power generated or consumed and (Idi is the flow 



constraint. Eq. (lei describes the physical properties of the 
network. In this formulation, p and v are the variables. Eqs. 
(Ici and (Id I are box constraints with respect to p which 



are the variables of the objective function (lai and they are 



relatively easy to handle. Eq. (Ibi together with (lei make 
the problem non-convex and hard to solve. To illustrate the 
algorithms, we first consider a simplified version of OPF with 
Ci2 — Cio — 0, Vi and neglect ( Ic i and ( Id i. Having — 
will not affect the optimal solution of the original problem. We 
will explain how to handle non-zero Ci2 later. By introducing 
a 71 X 71 complex matrix W — {Wik, 1 < i,k < n) ~ \\ 
we can write the simplified OPF in the sequel: 



H 



mmimize 



1=1 



CtlPi 



subject to 



rank(W^) = 1 

p = Re{diag(w^Y^)} 



(2a) 



(2b) 
(2c) 
(2d) 



Let C = diag(cii, C21, . . . , c„i) and M = (M^fe, 1 < i k < 
n) = i(Y^C + CY). By relaxing the rank constraint (29c 
we have the following semidefinite program (SDP): 



minimize Tr(MW) 
subject to 



w ^ 



(3a) 

(3b) 
(3c) 



where Tr(-) is the trace operator. We can solve this SDP 
at a central control center However, current algorithms for 
SDP, e.g. primal-dual interior-point methods |10|, can only 
handle problems with size up to several hundreds. We will 
decompose the problem into smaller ones by exploring the 
network structure. 

B. Zero Duality Gap 

By Q, the simplified OPF and SDP share the equivalent 
optimal solution provided that the network has a tree structure. 
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is a lossless cycle, or a combination of tree and cycle. For 
these kinds of network structures which are typically found 
in distribution networks, the optimal solution computed from 
(|3]l is exactly the same as that from d2|i. Targeting distribution 
networks, we can merely focus on (pjl. For completeness, the 
approach in ^4J is outlined below. 
The dual of ^ is given by 

n 

maximize ^{-XiV^^ + X^Vf) (4) 

i=l 

subject to Xi > 0, A, > Vz 

A + A/ ^ 0, 

where Aj and Aj are the Lagrangian multipliers associated with 
the constraints Wu < and y_^ < Wu respectively. From 
the KKT conditions, |4j showed that ^ always has a solution 
that is rank 1. 



C. Graph Structure 

We will use the following graph structures to decompose 
SDR 

Consider a graph G = {V,E), where V — {i\l < i < n} 
are vertices and E — {{i;k) £ V x V} are edges. Vertices 
i and k are adjacent if (i, k) E E. A clique C is a sub- 
set of V whose induced subgraph is fully connected, i.e., 
{i,k) G E,\/i,k G C. A clique is maximal if it cannot be 
extended to form a larger one by including any adjacent vertex 
to the clique. In other words, there does not exist a clique 
whose proper subset is a maximal clique. A chord is an edge 
which connects two non-adjacent vertices in a cycle. A graph 
is chordal if each of its cycles with four or more vertices 
contains a chord. Thus a chordal graph does not a cycle with 
four or more vertices. If G is not chordal, we can produce a 
corresponding chordal graph G = {V, E), where E ~ EU Ef 
and Ef = {{i,j) e V x V - E} are chords of G, called fill- 
in edges. G is not unique. From G, we can compute the set 
of all possible maximal cliques C — {Ci, . . . ,G\c\}, where 
Gi = {j e V} whose induced subgraph is complete and 
maximal. If G is a tree, each pair of vertices connecting by 
an edge forms a maximal clique. For a tree with n vertices, it 
can be decomposed into n — 1 maximal cliques. 



For M in ( 3a i, we can induce the corresponding G by having 
V = {i\l <i <n} and E ^ {{i,k)\M,^k ^ 0}. G has a very 
close relationship with the power network structure because of 
Y. If all c,i, Vi are non-zero, G directly represents the network. 

We use the following procedure to produce C from M. 

1) Construct a graph G = {V, E) from M. 

2) From G, compute Maximum Cardinality Search fTTl to 
construct an elimination ordering a of vertices |12|. 

3) With CT, perform Fill-In Computation p3] to obtain a 
chordal graph G. 

4) From G, determine the set of maximal cliques C by the 
Bron-Kerbosch algorithm p4) . 

Note that similar ideas about maximal cliques have been 
utilized to develop a parallel IPM for SDP y_5)-|[T7l. However, 
we make use of the ideas to decompose SDP into smaller 
problems, which can be tackled by any appropriate SDP 



algorithm, not necessarily IPM. Therefore, our approach is 
more flexible on that any future efficient SDP algorithm can 
be incorporated into our framework. 



III. Algorithms 

In this section, we will first present the primal and dual 
algorithms whose outputs are positive semidefinite matrices. 
Then we will explain the mechanism to convert such a matrix 
into the voltage vector. 



A. Primal Algorithm 

The objective function ([3a]l can be expressed as 

(5) 



Tr(MW) = MiW,,, 



i,k=l 

where each term M^^Wik can be classified into one of the 
following three categories: 

1) Ignored terms 

Each of which has M.^k = 0. Let I ^ {{i, k)\Mik = 0}. 

2) Unique terms 

For Mik 7^ 0, both i and k belong to a unique maximal 
clique. If e G/, then i,k ^ Gr,yr ^ I. Let U ~ 
{ii,k)\i,ke Gi,yi,i,k^ G^,Vr ^ I}. 

3) Shared terms 

For Mik 7^ 0, both i and k belongs to more than one 
maximal clique. 

Then (|5]l becomes 

Tr(MW)= M,",W,k+ E ^"W,^ 

i,k\{i,k)eX i,k\{i.k)eU-X 

+ Y MiW.k, (6) 

i,k\{i,k)^XUU 

where all the ignored terms can be ignored. Since each unique 
term is unique to each maximal clique, (|6]l becomes 



Tr(MW)= MiW.k+ Y 

i,fc|(i,fc)^XUW 



i,keCiyCiec 

\(i,k)eU-I 



(7) 



Eq. ( [3b] l gives bounds to each Wu, 1 < i < n and it is 
equivalent to 



<Wu<V,yieGi,VCi eC. 



(8) 



By [TSl, a matrix is positive semidefinite if all its subma- 
trices corresponding to the maximal cliques induced by the 
matrix are all positive semidefinite. Let Wc,c, be the partial 
matrix of W with rows and columns indexed according to G;. 



Eq. ( 3c I is equivalent to 



(9) 
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Hence, (|3]l is written as 

minimize ^ M^Wik 



i,keCiyCiec 
\(i,k)eu-x 



i,k\(i,k)fXUU 



subject to 



Wc,c, h oyCi e C. 



(10a) 

(10b) 
(10c) 



If we fix all Wik in the shared terms (those in the second 



summation in ( 10a i, ( 10 1 can be decomposed into \C\ subprob- 
lems, each of which corresponds to a maximal clique. For Ci, 
we have the subproblem I, as follows: 



rmmrmze 



subject to 



,keCi\{i,k)eu-x 



(11a) 



n <Wu<V,yieCi,iiCr,yr^l (lib) 
Wc,c, ^0. (11c) 

In ( [TT] ), only the semidefinite constraint ( |1 lc| i involves those 
variables which are not unique to Ci, i.e., Wik such that 
(i, k) ^ U. By introducing a slack variable Xi^j = Wik 
for each shared Wik in subproblem I, we define Wc,c, = 
(l^jfe, i, fc e C;) and Mqc, = (M^fc, i,keCi) where 



and 



W,k 
Xik,i 



M,k 




if {i,k) e U 
otherwise 



if [i,k) eU 
otherwise. 



Then ([TT]) becomes 

minimize Tr(Mc,CzWc,c,) 
subject to 

a2 



(12a) 



Vf <Wu<V,yieCi,ii Cr^r ^ I (12b) 
Wc,c, ^ (12c) 
X,k,i^W,k,'^i,k\{i,k)iU. (12d) 



Note that Wife's in (12dl are given to the subproblem. 
When given such Wife's, all subproblems are independent 
and can be solved in parallel. Let the domain of ([12]) be 
Given W,k where {i,k) ^ U, let 4)i{W,k\{i,k) i U) = 



inf. 



._{Tr(Mc,c,Wc,c, 



becomes 



minimize ^ (l)i{W,k\{i,k) ^U) + E ^^i^^^fc 

VCiGC i,k\(i,k)<$_U 

(13a) 

subject to 



(13b) 



Eq. ( [T3] l is the master problem which minimizes those Wik 
shared by the maximal cliques. With those shared Wik com- 
puted in ( [T3] l, we minimize the Wik unique to each subproblem 
given in (|12[). 



In ( [T3] l, those shared Wife can be further classified according 
to nodes and edges: 
1) Nodes 

Let Xii i be the Lagrangian multiplier for (12di with 
i ^ k. The subgradient of Wu with respect to sub- 
problem I is —Xii^i 1191. Thus the overall subgradient is 
I]/Kec, + iteration t, we update Wi^ 

by 



Wl-+^' = Proj I W^V' - " 




where 



Proj(x) 



n 


if a; < V^, 




if a; > 




otherwise, 



a*^*^ is the step size at iteration t, and W^^' represents 
Wii at iteration t. 
2) Edges in E 

We consider (12di with i ^ k. Since Wik and Xife j 



are complex numbers, we can handle the real and 
imaginary parts separately, i.e. Re{Xife.i} = RejWife} 
and Im{X,fe,j} = Im{W,fe}. Let Af^= ; and Aj^; be their 
corresponding Lagrangian multipliers for subproblem /, 



respectively. In ( 20a i, Vi 7^ k, the ik and ki terms always 
come in a pair. We have 

MgW^k + Mfe^Wfe, =2Re{M^}Re{W,fe} 

- 2Im{M,fe}Im{W,fe}. 

A subgradient of the real part of Wife is 
E/K,feec, (-Affe^,;) + 2Re{M^}. At iteration t. 



we update Re{W-'^^} by 
Re{W,f^)} =Re{W,<*)}- 



v(*) 



E (-Affe^,)+2Re{M^} 

\l\i,kl^Ci 



(15) 



Similarly, for the imaginary part, we have 

lm{wll+'^} =Im{W,<*)}- 



v(*) 



E (-AflJ - 2Im{Mf } 
a\i,k^Ci 



(16) 

3) Edges in Ef 

Recall that fill-in edges are "artificial" edges added to G 
to make G. For (z, k) e Ef, we have Mik = 0,i ^ k. 
Similarly, at iteration k, we update its real and imaginary 
parts by 



Ro{wt'^} = Rc{W^l^} 




(~Ajfe ,) 



(17) 
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and 



Im{M^«} 




(18) 



We can interpret the updating mechanism as follows: certain 
maximal cliques share a component Wit (if i ~ k, it 
corresponds to a node; otherwise, it corresponds to an edge 
or a fill-in edge). Wik represents electricity resources and 
—M^ is its default price. An agent (i.e. a node responsible for 
computing the update) which is common to all those maximal 
cliques sharing the resource determines how much resource 
should be allocated to each maximal clique. In other words, it 
fixes Wik and every party gets this amount. Xik^i is the actual 
resource required by Ci and Xik^i corresponds to the price 
of the resources when Wik is allocated to it. If C; requires 
more resource than those allocated, i.e., Xi^ i > Wik, then 
Xik,i > 0. If the net price, i.e. J2i\i.kec, ^ik.i-Mik, is positive, 
the agent should increase the amount of resource allocating to 
the maximal cliques because it can earn more. If the net price 
is negative, then supply is larger than demand and it should 
reduce the amount of allocated resources. 

From (14i-(18i, all shared Wik can be updated indepen- 



dently. The update of each Wik only involves those {C;|C; G 



C,i,k e Ci}. In other words, (13i can be further computed 
separately according to those maximal cliques shared by each 

W,k. 

The pseudocode of the primal algorithm is as follows: 

Algorithm 1 Primal Algorithm 

Given Q,V,VX 
1. 
2. 



12 1 



Construct (|12| for each maximal clique 
while stoppmg criteria not matched do 
3. for each subproblem I (in parallel) do 

4. Given Wik with {i,k) ^ U, solve 

5. Return Xik^yi, k\{i, k) 

6. end for 

7. Given A,;j. ;VZ|i, k G Ci, update the shared Wik with 
(14])-([l8l (in parallel) 

enawhue 



B. Dual Algorithm 

Let Vlik = {Ci\i,k e C/,V?}. Problem ^ can be written 

as 



Tr(MW) 



E 

iM{i,k)eu 



MgWik 



J2 

i,k\{i,k)^U 



\^ik\ 



E 

c,ec \i,keCi\{i,k)eu 



E ^^"w.>^ 



E 

t,k£Ci\{i,k)(f_U 



MgW^k 



'ik\ 



(19) 



Problem ([3]) becomes 
/ 



mmimize 



E 

. i.keCi 



M"kW,k 



MgWik 



\ 



i,keCi 
\{i,k)fU 



I a 



subject to 



<Wu<vlyieCiyQeC 



w, 



CiCi 



h 0,VC, eC. 



/ 

(20a) 

(20b) 
(20c) 



Problem (20i can be separated into subproblems based 
on the maximal cliques. However, the subproblems are not 



completely independent of each other as ( 20c i involves some 
common variables shared between the subproblems. Similar to 
the primal algorithm, we can replace WciCi with W^Cr For 
each Wife I (i, k) ^ U, let Xik^i be a copy of Wik in C; £ ftik- 
To make all Wc,Ci consistent, we should have 

Xi, 



W,k = X 



ikAi 



^ik,l2 



X. 



ikJii 



yw,k\{i.k)iu, 

(21) 

or simply 

Xik^li — Xik,l2 = ■■■ = -''^ifc./is-i.^i ,V^r|Cir G fiife, 

yi,k\{i,k) (^U. (22) 

For each {i,k) ^ U, (22 1 can be written into \ilik\ — 1 
equalities, e.g.. 



XikJi — Xik.l2 
Xik,l2 — Xik.U 



(23) 



X, 



k.l. 



X, 



k.lu 



As shown later, the update mechanism of the dual algorithm 



depends only on how we arrange (22i into equalities. In fact, 
there are many ways to express the \Q.ik\ — 1 equalities pro- 
vided that each Xik^i appears in at least one of the equalities. 
Suppose the rth equality be Xik,r{^) ~ X^/j r(2). We assign 
a Lagrangian multiplier Vik,r to it. Then we have 



Vik,i{Xiks{^) 

Vik.2{,Xik^2{^) 



x^kAi)) 

X,fc,2(2)) 



(24) 



'Uik.\ai^\~i{Xik,\Q.ik\-iW - Xik.\aik\-i{'^)) = 



When we sum all these equalities up, each Xik^i will be as- 
sociated with an aggregate Lagrangian multiplier Vik.i, which 
is composed of all Vik associated with Xik.i- For example, 
in (|23]l, we have Xiksi'^) = X,kM^ Xik^i{2) = X^kM^ 
and Xik{2,l) = Xik.i^- Thus Vik,u = ""ik^i and v^kM = 

Vik,2 — Vik,l- 

In ( [24) i, yi < i, k < n, the corresponding rth equality for 
the {i,k) pair implies the one for the {k,i) pair, i.e., 

Xik,r{^) = Xik^r{2) Xki^ri^) = Xki,r{'2), 
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due the positive semidefinite property of W given in ( [3q l. We 
have 

=^ (Vik^rXik^ri^))^ = (Ui/c,r-^ifc,r(2))^ 
^v"k,rXg,rW=vg^,Xg^,{2) 

Thus the aggregate Lagrangian multipher for X^i^i can be 
computed directly from that for Xik.i, i.e., Vki^i = vg 

Let V = {v,kd,.,i,k £ Ci\{i,k) ^ UyCi £ C;lr\Ci^ £ 
ilik)- We form the dual function d{v,W) by aggregating (22 1 
into ([T9|. We have 

d(u,W) 



Let v^*/}^, a^*) > 0, and x|^.\(z) be the Lagrangian multiplier 
of the rth equality associated with Wik, the step size, 
respectively, at time t. Then we can update Vik,r in (p4)l by 



,(*+!) 



If ^Ik.ri^) < Kk'.ri^l then v.^ 



,(*) 



(28) 



(t) 



> 4k.r- This will 
make the coefficient corresponding to Xik^ri^) larger while 
making that corresponding to Xik^r{2) smaller. At time t + 
1, the subproblem will obtain xllX^\l) < xll^^l) and 



xlZ^^m > Xll'm. Hence, \XZ^'(2 



C^(t+i) 



ik,r \ 



xr.'>{i)\ < 



E E 

Ciec \i,k£Ci\{i,k)£U 

E E ^ik,l^Xik,l^ 



E 

i,k£Ci\{i,k)(^U 



ik 



IXll'A^) - ^ikM)\- 0° the other hand, if X^l>A2) > 
■^ik,ri^)' then Ujfc^^^ < t^-fc^r- This will make the coefficient 
corresponding to Xik,r{^) smaller while making that corre- 
sponding to Xik^r{'2) larger. Then we will get X^l^J'\l) > 



ik\ 



Xl^m and Xir'r>{2) < XfAi"^). This wiU also make 



■r{t + l) 



it) 



\x. 



(t + 1) 



(2) 



{2)~XIII{1)\. Therefore, 



(28 1 drives Xik/s in (|22|) become closer to each other in value 



when the algorithm evolves. In other words, (28 1 tries to make 



E 

Ciec 



1 i,keCi 

\\(t:k)eu 



ik 



E 

i,keCi 

\(i,k)^U 



ik 
\^ik\ 



equality (22 1 hold when the algorithm converges. 
\ At any time before the algorithm converges, i.e.. 



Ciec 

Given v, (|20ll becomes 



^ikjXik.l 



(25) 



(22 1 does not hold, the solution W with the computed 
Xife^;, Vi, fc) ^ Uyi\Ci £ ilik is an infeasible solution. 
/ We can always construct a feasible W with Wik which is the 
average of all Xn^ j in (22i. 

The purpose of (28 1 is to make the two entity -^^^''^(l) 



and Xl'^^{2) closer to each other As long as ^^^.''^(l) and 
xff}^{2) have been computed (from two subproblems), we 



minimize (i({;,Wc,Ci) 



CiSC 



subject to 



(26a) 



(26b) 
(26c) 



can update Vik^r with 28 Thus different Vik can be updated 



Then problem ( 26 1 can be divided into subproblems according 
to the maximal cliques and each of them is independent of 
each other Subproblem I is stated as: 



minimize ^ MgW,k+ ^ 



i,k£Ci 
\{i,k)eU 



i,keCi 
\{i,k)<^U 



^" +7, \X 
77;, — r + Vikj I Aik,i 

"ife 



(27a) 



subject to 

< <v\Ai(^Ci\{i,i)£U, (27b) 

Y}<Xuj<vlyi£Ci\{i,i)iU, (27c) 

Wc,c, > 0. (27d) 

By solvin g (|27| , we denote the optimal Xn.,r{z) for the rth 
equality in (^24jrby X°^^^{z), where z £ {1,2}. The gradient 
of — W(7,Cif]with respect to Vik^r is 

-^=XrA2)-XZ{l)- 

^^ik,r 

'in the dual form, we maximize inf^ii(i5,W) over v. In minimization, 
we consider ~d{v,V/c^c^). 



asynchronously. Since the only co-ordination between sub- 
problems is through ( [28] l, synchronization is not required in 
dual algorithm. 

We can interpret the updating mechanism as follows: Vik.r 
is the price assigned to equality Xik^r{^) = Xik^r{2). We 
can treat Xik^ri^) and Xik.r{2) as demand and supply of 
electricity resources, respectively. If the demand is larger than 
the supply, i.e., X^fe ,.(1) > Xik,r{2), we should increase the 
price so as to suppress the demand and to equalize the supply 
and demand. On the other hand, if the supply is larger than 
the demand, i.e. Xik.r{i) < Xik^r{2), we should reduce the 
price in order to boost the demand. 

The pseudocode of the dual algorithm is as follows: 

C. Computation of Voltage 

When either the primal or the dual algorithm converges, 
assuming zero duality gap, we obtain the optimal W = vv^. 
To obtain each bus voltage and voltage flown on each line, 
we first compute the voltage magnitude at each bus, \Vi\ = 
^/Wii, 1 < i < n. For 1 < i, k < n, if there is a line between 
nodes i and fc, the corresponding line voltage angle difference 
9ik can be found by solving Wik = \Vi\\Vk\e^^^'' at either 
bus i or k. For the former, bus k needs to send \Vk\ to bus 
i, and vice versa. By fixing the voltage angle of a particular 
bus to zero, the voltages of the whole network can be found 
subsequently. 
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Algorithm 2 Dual Algorithm 



Given Q,V,V,C 

1. Pair up slack variables for the shared variables into equalities 

2. Construct ( 27 1 for each maximal clique 



while stoppmg criteria not matched do 

4. for each subproblem I (in parallel) do 

5. Given Vik, solve (27i 

6. Return Xik^i^i, fc[(i, k) <^U 
1. end for 

5. Given Xik.i^, update the price Vik,i^ with (28 i (in parallel 

and asynchronously) ^'g- ^■ 

end while 




(a) Toplogy 



n-bus radial network 



Primal 



Dual 



(b) Communications between 
nodes I and n 



IV. Quadratic Cost Function 

Up to now we have focused on the OPF problem with a 
linear objective function. In practice, sometimes a quadratic 
cost function is used. If this is the case, the methods developed 
so far can be used as subroutines to solve the OPF problem 
by adding a outer loop to the iteration. 

Let costi(Pi) = Ci2P^ + CiiPi be the cost function as- 
sociated with Pi. We assume this function is convex for all 
buses, that is, Ci2 > Vi. From ( fTe] ), Pi ~ Tr(AiVv^), where 
Ai = l{{Y"Ei) + EiY) and Ei is the matrix with 1 in the 
{i, i)th entry and zero everywhere else. Now the OPF problem 
is (compare with (|2]i) 



minimize Ci2Tr(A,VF)^ + caTr(AiW^) (29a) 
subject to 



rank(M^) = 1 

p = Re{diag(vv^Y^)} 



(29b) 
(29c) 
(29d) 



Using Schur's complement, we may write p9l) equivalently 



as 



minimize ^y^^jtj + CiiTr{AiW) 



subject to 



U Vc^Tr(AiWy 
/^Tr(AiW) 1 

< Wu<Vi\yt 



(30a) 



^ Vi (30b) 



rank(W^) = 1. 

Relax the the rank 1 constraint and taking the dual, we get 

n 

^2 



maximize AiVj + A^F^ — Ui) 

1=1 

n 

subject to ^^(cii^j — '2^/a^ZiAi) + A ^ 



'i:=l 
1 

Zi Hi 



OVi, 



(31a) 



(31b) 



(31c) 



where the constraint pic[ ) corresponds to the Schur's compli- 
ment constraint in (|30]l. The constraint pic| i can be rewritten 
as Ui > zf, for a given Zi, the maximizing Ui is zf. Therefore 



the we may drop the constraints ( |31c| l and replace the Ui in the 
objective function by zf. If we fix the z/s, then ( |3T| i becomes 
a function of z. 's 



J(z) = maximize ^(-A,;V^- + X,V^ - zf) (32) 

i=l 
n 

subject to ^(cii^i — 2y/ci2ZiAi) + A ^ 0. 



For fixed z, ( (32[ l is in the form of Q. Therefor J(z) is a dual of 
the optimization problem with linear cost functions with costs 
(cii - 2yc]Izi,C2i - 2^/c2^Z2, . . . ,c„i - 2^/c;;2^z„). To find 
the optimal solution of ^32\ we may use any of the algorithm in 
the previous sections. Let W*{z) denote the optimal solution 
to J(z). To find the optimal z, we use a gradient algorithm. 
The Lagrangian of ((32]i is 



£(A,W^)=5;](-A,F- +A,;Z?-^?) 

n 

+ Tr((^(QiA:-2y^z,A,)+A)W^). (33) 



By a standard result in convex programming, the gradient of 
J{z) is given by ^ = "^^j^ = -JYxi^^AS^*) ~ 
2zi, where {\* ^W*) is a pair of optimal dual-primal solu- 
tions (dependent on z). Therefore, to solve the problem with 
quadratic cost functions, we add an additional outer loop to 
the solution algorithms for the linear cost functions. 



V. Illustrative Examples 

In this section, we will consider two examples to illustrate 
how the algorithms work. The first one is a n-bus radial 
network with height equal to one while the other is a five- 
bus network with a four-bus ring. The former gives ideas 
how many maximal cliques share a common bus. The latter 
demonstrates multiple-level bus and edge sharing. We will also 
show where the different components are implemented and 
how the communication is accomplished. 
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A. n-Bus Radial Network 



Consider the topology of the network given in Fig. 1(a) We 
have 



(Mil 













M22 













M33 • 


■ Man 



M = 

\Mni M„2 M„3 ••• M„„y 
1} Primal Algorithm: Eq. (j7]l becomes 



Tr(MW) 



n-l 

1=1 



shared term 



unique terms 



Each branch with the end nodes forms a maximal clique. We 
have C = {C/|l < I < n - 1} where Q = and U = 

{n,n). By introducing a slack variable X„„./ for Wnn to C;, 
we have 



Let 



Wnl Xnn.l 



( Mil Ml, 



\Mnl 

Given W^„„, subproblem I for C; is stated as 

Tr(M(Wc,c,) 



mmimize 



subject to y_i < Tr 



1 




W, 



c,c. 



< V, 



Trig JjWc,c, =W^nn 
Wc,c, h 



which is an SDP with a 2 x 2 vaiiable. Recall that (?!)/(W„„) 
is the optimal value of subproblem I given Wnn- The master 
problem is 

minimize J2 ?=i 0i ( W^nn ) + M^ Wnn 
subject to y2 < T^^^ < 

A„„_; is the Lagrangian multiplier for equality Xnn,i = Wnn 
of subproblem I. We update Wnn by 

Wii^'^ = Proj (^W^W - (^g (-A„„,) + Mf„ j j 

(35) 



In iteration bus n broadcasts Wnn to bus 1,1 < I < n — 1 
After receiving W,l5], For all I, bus ^ solves its own subproblem 
(|34]l by any suitable SDP method, e.g. primal-dual IPM |10|, 
in parallel and then returns A„„ ; to bus After receiving 
all Xnn,i, node n updates Wnn by (35 1. The communication 



pattern is shown in Fig. 1(b) 



^For most of the interior-point metliods, primal and dual solutions come in 
pair. When the algorithm finds the optimal Wo, Ci> it will also give A„„ j. 
Thus no extra calculation is required to determine \„„ 



2) Dual Algorithm: Eq. (fT9]l becomes 



TriMW) = J2 (MnlWnl + M^^Win + MllWii 



1=1 



M^LWnn. 

n — 1 



As only Wnn is common to all maximal cliques, we have 
^^nn = {Ci, . . . , Cn-i} and 



(36) 



Assume that ( 36 1 is arranged as follows. We assign Lagrangian 
mulipliers to the equalities and we have 

1Jnn,liXnn,l — Xnn,2) — 
Vnn,2(Xnn,l ~ Xnn,?,) — 



(37) 



For 1 < ; < rt - 2, X„„_;(l) = Xnn,i and X„„,;(2) = 
Xnn,i+i- Then we have 

^nn.l — '^nn,l ^" ' ' ' ^" ^nn.n— 2i 
W„n./ = -WnnJ-l, 2 < / < n - 1. 



Let 



where 



M, 



Mil Min 
Mnl Ml 



(34) Mi^ 




+ Vnn,n-2j 7 ^ — 1 



2 < ? < n - 1 



and 



Subproblem Z is stated as 
minimize 



Wll Win 
Wnl Xnn,l 



Tr(M,Wc,c,) 



subject to Z?<Tr(J [] ) Wc,c, 



< Tr f [j J ) Wc,c< < Vl 



(38) 



which is also an SDP with a 2 x 2 variable. We update the 
price by, for 1 < r < n — 2, 



,(t+i) 



a''*'' (-'^nn,r+l " Xnn.l) 



(39) 



At time t, node Z from Cj,l < / < n — 1, solves (38|p 



e.g., by IPM, in parallel and sends X„„, from the optimal 



WciCi to bus n. Whenever any pair of Xnn specified in (47i 



(e.g., Xnn,i and Xnn,i for Ci and C;, respectively) reach bus 
n, price Vnn.i-i can be updated with (39 1 by bus n and the 
updated Vnn,i-i is multicast back to the corresponding buses 
(e.g., nodes 1 and 0- The communication pattern is shown in 



Fig. 1(b) 
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(a) Primal algorithm 




«44,1 

(b) Dual algorithm 

Fig. 2. Structure and communication patterns for the five-bus network with 
a four-bus ring 



B. Five-Bus Network with a Four-Bus Ring 



Consider the topology of the network given in Fig. 2(a) We 
have 



M 



To make it chordal, suppose we add a fill-in edge between 
buses 2 and 3 and we have Ci = {1,2,3}, C2 = {2,3,4}, 
and C3 — {4, 5}. Assume bus 2 is used to co-ordinate Ci and 
C2, and bus 4 for C2 and C3. 

1) Primal Algorithm: Let 





M12 


Mi3 





^ 


M21 


M22 





M24 





M31 





M33 


M34 








M42 


M43 


M44 


M45 


I 








M54 


M55J 



'Mil 
M21 






A/12 Ml3^ 





M43 



M54 





Af24^ 
Af34 





WciCi = I ^^21 ^22.1 ^23.1 1 , Ml 

WC2C2 = I ^32,2 ^33,2 W3i ) , M2 
y W42 H^43 -''^44,2 / 

WC3C3 - w,, W,s) ' 
We have 



'in fact, any bus in a maximal clique can be elected to solve the subproblem. 
In this case, node I is chosen to reduce the computation concentrated at bus 



A/45 
A/55 



Subproblem 1: Given W22 and W33, 

minimize Tr(MiWciCi) 
subject to Vl <Wii < v\ 

WciCi > 0. 
Subproblem 2: Given W22, ^33, and W44, 



(40) 



minimize 



subject to Xi, 



Tr(M2W, 

.2 = Wu 



C2C2) 



2,3,4 



(41) 



W, 



O2 C2 



>- 0. 



Subproblem 3: Given W44, 



minimize Tr(M3Wc'3C3) 
subject to Vl < W55 < V\ 

^44^3 — W44 

WC3C3 > 0. 



(42) 



In iteration t, bus 2 sends W22 and W23 to both buses 1 and 
4 and bus 3 sends 1^33 to both buses 1 and 4. Bus 4 sends 
W44 to bus 5. Buses 1, 4, and 5 compute (|40]), (|4T}, and ( [42] i, 
respectively. Then bus 1 sends A22,i to bus 2 and A33.1 to bus 
3. Bus 4 sends A22.2 to bus 2 and A33 2 to bus 3. Bus 5 sends 
A44,3 to bus 4, which has A44_2- Buses 2 and 3 update Wa by 



(t+i) _ 



Proj ( 



it) 



M„ 



(43) 



within [V_^, V.^], for i = 2, 3, respectively. Bus 4 updates VF44 
by 



(t+i) _ 



Proj ( 



it) 



v(*) 



( — A44^2 — A44^3 + M£^^ 



within [Vl, Vl]. Bus 2 updates W23 by 



Re{M^; 



(t+i)i 



23 } = Re{VK2^^a - a 



(-A 



Re 
23,1 



(44) 



(45) 



and 



Im{T^i*+^'} = MT^i*)} - (-A'™,i - X'^,^2)- (46) 



The communication pattern is shown in Fig. 2(a)| 

2) Dual Algorithm: Assume the equalities for the slack 
variables are arranged as follows. With Lagrangian multipliers, 
we have 



(47) 



f 22,1 (-^22,1 


— X22.2) 


= 


1^23,1 (-'^23,1 


^ -^^23, 2) 


= 


^'33,1 (-'^33,1 


— -^^33, 2) 


= 


f44,l(-''^44,2 


— ^44^3) 


= 



'Let 




+ ""22,1 
,H 
23,1 



Mi3 

^ +^^23,l 

4^+^^33,1 



^^-^^23.1 
^-t^33.1 
M43 



M24 
M34 
^ + ^^44,ly 
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TABLE I 

Normalized CPU time for distribution test feeders" 



Number 
of buses 


Centralized 


Cumulative 


Distributed 


Primal 


Dual 


Primal 


Dual 


8 


1.85 


7.21 


5.62 


1.52 


1.00 


34 


298.68 


37.79 


33.89 


1.94 


1.70 


123 


_b 


143.39 


126.48 


2.24 


1.64 



^ The CPU times are normalized by 0.0857s. 
The solver carmot be applied because of the out-of~memory problem. 



and 



M54 



M45 
M55 



We have 

• Subproblem 1: Given V22,i, 1^23,1, and ^33,1, 
minimize Tr(MiWciCi) 

1,2,3 



subject to <W.u<V^i,i 



(48) 



W, 



CiCi 



y 0. 



Subproblem 2: Given ^22,1, "^23,1, ^33,1, and ^44,1, 

minimize Tr(M2 C2 ) 

subject to V^^ <Wu<V^i,i^2,3,A 
Wc.c. h 0. 
Subproblem 3: Given U44 1, 

minimize Tr(M3Wc3C3) 
subject to V] < Wu <v],i = 4,,b 
WC3C3 > 0. 



(49) 



(50) 



At time t, bus 2 announces U22.1, "^23,1, and U33.1 to buses 
1 and 4, and bus 4 sends U44.1 to bus 5. Buses 1, 4, and 



5 solve (48 I, (49 1, and (50 1, respectively. Then bus 1 sends 
-^22,1, ^23,1, and X33 1 to bus 2. Bus 4 sends X22,2, -''^23, 2, 
and X33 2 to bus 2. Bus 5 sends X44.3 to bus 4. Bus 2 updates 

W22,i' "^23,1 and W33,i by 



(t+1) _ (t) 



X^k,i) ,i,k = 2,3,« < fc. 



and bus 4 updates W44.1 by 



,(*+!) 



,(*) 



(-''^44, 3 ~ -^^44, 2) ■ 



(51) 



(52) 



The communication pattern is shown in Fig. 2(b)[ 



VI. Simulation Results 

To evaluate the performance of the algorithms, we perform 
extensive simulations on various network settings. Since OPF 
is formulated as an SDP, the optimal solution can be computed 
in polynomial time by any popular SDP algorithms, e.g. IPM. 
Recall that the primal or dual algorithm aims to divide the 
original problem into smaller ones and to coordinate the 
subproblems, which of each can be solved by any SDP solver 
independently. The primal and dual algorithms perform coor- 
dination by indicating what problem data should be allocated 
to each subproblem and do simple calculations to update the 
shared terms (for the primal) and the prices (for the dual). 



TABLE II 

Normalized CPU time for IEEE Power transmission system 

BENCHMARKS" 



Number 


Centralized 


Cumulative 


Distributed 


iterations 


Initial 


of buses 




dual 


dual 




step size 


14 


5.38 


5.38 


1.00 


1 


30 


30 


45.29 


58.60 


5.38 


6 


30 


57 


1696.79 


49.08 


4.28 


4 


30 


118 


_b 


704.46 


13.51 


9 


300 



The simulations for this problem set are done on MacBook Pro with 2.4 GHz Intel 
core 15 and 4 GB RAM. The CPU times are normalized by 0.1410s. 
The solver cannot be applied because of the out-of-memory problem. 





-centralized 


■ ■ ■ o- 


■ cumulative primal 


■ ■ ■ G- 


■ cumulative dual 


— -1- 


-distributed primal 


— *- 


-distributed dual 



^ 10 

Ill 
E 

10" 




Fig. 3. CPU time of the various approaches on radial networks with bounded 
voltages 




20 30 40 

Number of buses 



50 



Fig. 4. Success rates of the primal and dual algorithms on radial networks 
with bounded voltages 



When compared with those done by the SDP solver, the 
computation and ordination required solely by our algorithms 
are relatively far less stringent. As a whole, the bottleneck of 
computation should be at the SDP solver In our simulation, 
we program the primal and dual algorithms in MATLAB and 
and solve each SDP with YALMIP |20| and SeDuMi |21|. 
To get rid of the dependence on the programming language 
and to simplify the comparison, we only count the CPU 
time spent on the SDP solver. Moreover, we can arrange 
the subproblems to be solved in a single node or distribute 
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them to different nodes in the network. For the former, we 
assume the problems are handled sequentially and we call it 
the cumulative approach. The latter, named as the distributed 
approach, addresses the subproblems in parallel. Without our 
algorithms, the (original) problem will be solved in its original 
form in a centralized manner. Here we compare the CPU 
times required for the SDP(s) among the centralized approach, 
(primal and dual) cumulative approaches, and (primal and 
dual) distributed approaches. 

We run the simulations on Dell PowerEdge 2650 with 
2 X 3.06GHz Xeon and 6GB RAM (except those for the 
transmission system benchmarks in Table In order to 
monitor the performance in each simulation run, we assem- 
ble the partial solutions (done by the subproblems) to form 
a complete one for the original problem and evaluate the 
corresponding objective valuej^ Our algorithms stop when 
the computed objective value falls in the range of 10^^ x 
the global minimum. We assume that the dual algorithm is 
synchronized. In other words, all subproblems for the dual 
are solved in each iteration (but this is not required when 
implemented in real systems). The initial step size a'"' is set 
to one and it is updated by a^*^ = a'^*-^'^^ /t^\/t > 0. An 
algorithm is deemed successful if the stopping criterion is met 
in 100 iterations. 

We perform simulations on three problem sets; the first two 
focuses on tree-like networks while the last one is about trans- 
mission networks. The first problem set is some distribution 
test feeder benchmarks p2) . As the data set does not specify 
the cost function of power production/consumption, we create 
a problem instance by randomly generating the costs. To do 
this, we first select one node, e.g. node i, to be the power 
source node with cn set randomly in the range (0, 10). For 
other node k ^ i, Cki set randomly in the range (—10,0). 
We create 100 instances for each network. Table U shows 
the averages of the normalized CPU times of the various 
approaches. All algorithms converge in 100 iterations for all 
instances. 

The second problem set is the n-bus radial network demon- 
strated in Section [V] For each instance, the root is the power 
source with a random cost selected in (0, 10) and each of 
the rest takes a random cost in (—10,0). For each node i, 
we specify a number ^ in (0.9, 1.1) and set — 0.95^ and 
Vi = 1.05^. For each line, the magnitudes of the conductance 
and susceptance are randomly assigned in (0, 10). We produce 
100 instances for each n and plot the average CPU times in 
logarithmic scale in Fig. [3] 

From the simulation results for these two problem sets, both 
the primal and dual algorithms converge very fast for distri- 
bution networks. The CPU time for the centralized approach 
grows very fast with the size of the network. The CPU time 
grows roughly linearly for the cumulative approach while it 
becomes independent of the size for the distributed approach. 
We define success rate as the fraction of the total number of 
simulation runs with stopping criterion met in 100 iterations, 
shown in Fig. |4] The success rate of the primal is almost 

The results in Tables |l]and|ll] are normalized, and thus, they are compa- 
rable. 

^The assembly of partial solutions is not required in real implementation. 



100% for all tested network sizes. The dual fails to converge 
in 100 iterations for a small fraction of small networks but 
the success rate grows to almost 100% with the network 
size. In general,the primal and dual algorithms are similar in 
performance but the dual requires a little bit less CPU time 
than the primal on the average. 

The third problem set is some IEEE power transmission 
system benchmarks f23l. As pointed out in f3l, these test 
cases have zero duality zap although they have network 
structures different from what we mention in Section III-BI 
Table |ll] shows the CPU times required, the iterations for 
convergence, and the initial step sizes. The primal algorithm 
is not applied to this problem set and the reason will be 
given in the next section. For these transmission network 
topologies, the maximal cliques of the fill-in graphs are much 
irregular than those with tree-structured networks. There are 
many ways to construct the maximal cliques and different 
construction can result in different convergence speed. The 
study of the relationship between maximal clique construction 
and the algorithm performance is out of the scope of this paper 
In this simulation, we randomly choose one maximal clique 
configuration and the step sizes are adjusted individually so 
as to have fast convergence. Nevertheless, the dual algorithm 
is more desirable than the centralized approach. 

VII. Discussion 

Both the primal and dual algorithms try to tackle the original 
problem by solving smaller subproblems but the ways to 
handle the information corresponding to the common partial 
solutions between subproblems are different. For the primal 
algorithm, if a common solution corresponds to a bus, that bus 
will compute its partial solution with the required information. 
For example, in Fig. 1(a) Wnn for bus n is common to all 
the subproblems. Bus n computes Wnn with its own y„, Vn, 
and Mnn- Only the computed Wnn is required to transmit to 
other buses which do not require bus n's information. For the 
dual algorithm, each subproblem needs to acquire all its bus 
information, even for the common bus. Consider the example 
in Fig. 1(a) again, for 1 < i < n — 1, node i requires 
V_n, Vn, and Mnn to solve its subproblem. Similarly, if the 
common solution is for a link, we can just elect any one of 
the connected bus to do the computation with the primal and 
dual algorithms. Therefore, the primal algorithm requires less 
information sharing between buses and it favors situations with 
sensitive bus information. 

Our algorithms do not require an overlay of communica- 
tion networks with different topology. From the examples in 
Section |Vj all the communication is one-hop. A node needs 
to transmit data, e.g. W, A, X, and v, to its neighbor nodes 
only. Communication links need to be built along with existing 
transmission lines only. 

The primal algorithm is suitable for networks with tree 
structure while the dual can handle those with rings. In fact, the 
primal is not very efficient to update the partial solutions for 
(fill-in) edges, i.e. ([T5]l-( 18 1, especially when their values are 
closed to the boundary of the feasible region. When updating 
the variables, the bus variables, i.e. Wu, are bounded but it is 
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not the case for the edge variables, i.e. Wik,i 7^ k. When the 
step size is too large, (15i-(18i may drive the current point 



out of the feasible region. If so, the obtained A cannot be 
used to form a subgradient, and thus ( 15 i-lfTSjl 



fail. When 

reaching an infeasible point from a feasible one, we can always 
step back and choose a smaller step for an update again. 
When approaching the optimal solution which is located on 
the boundary for SDP, this step-back procedure is not very 
efficient. The dual algorithm does not have this problem when 



updating v with ( 28 1. As mentioned, we can always constrain 



a feasible solution by averaging the shared variables. 

Tree networks have a very nice property with respect to 
maximal cUque construction; each branch with its attached 
nodes form a maximal clique. Hence, it is straightforward to 
decompose the problem into subproblems. However, when it 
comes to a more irregular network like those IEEE transmis- 
sion system benchmarks, the number of ways to decompose 
the problem grows with the size of the network. The perfor- 
mance of our algorithms also depends on how we form the 
maximal cliques and the initial step size needs to be adjusted 
accordingly. For problems with tree structure, the performance 
is easier to predict and the algorithms converge faster. 

The dual algorithm is more resistant to communication 
delay than the primal. For the primal, an update of a shared 
variable requires A from all involved subproblems and thus 
synchronization is required. Delay of computing or transmit- 
ting an A from any subproblem can affect the whole algorithm 
proceed. On the other hand, an update of an v requires the 
X from two pre-associated subproblems according to the 



arrangement of the inequalities in 22 Delay of computing 



or transmitting a particular X can affect the update of some 
but not all the v. Thus the dual algorithm is asynchronous. 



Moreover, we can pair the variables in ( 22 1 into equalities 
differently and secretly whenever we start the dual algorithm. 
In some sense, the dual algorithm is more robust to attack 
stemmed from communication on the communication links. 



VIII. Conclusion 

OPF is very important in planning the schedule of power 
generation. In the smart grid paradigm, more renewable energy 
sources will be incorporated into the system, especially in 
distribution networks, and the problem size will also grow 
tremendously. As problems with some special structures (e.g. 
trees for distribution networks) have a zero duality gap, we can 
find the optimal solution by solving the convex dual problem. 
In this paper, we propose the primal and dual algorithms (with 
respect to the primal and dual decomposition techniques) to 
speed up the computation of the convexified OPF problem. 
The problem is decomposed into smaller subproblems, each 
of which can be solved independently and effectively. The 
primal algorithm coordinates the subproblems by controlling 
the shared terms (related to electricity resources) while the 
dual one manages them by updating the prices. From the 
simulation results for tree-structure problems, the computation 
time grows linearly with the problem size if we solve the 
decomposed problem in a central node with our algorithms. 
The computation time becomes independent of the problem 



size when the subproblems are solved in parallel in different 
nodes. Even without nice network structure such as a tree, the 
dual algorithm outperforms the centralized approach without 
decomposition. Therefore, the primal and dual algorithms are 
excellent in addressing OPF, especially for distribution net- 
works. In future, we will improve the algorithm by incorporate 
more constraints into the OPF problem and move to nonlinear 
objective functions. 
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